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Abstract 

■ The presence of a chemical potential completely changes the analytical structure of 

q . the QCD partition function. In particular, the eigenvalues of the Dirac operator are 

distributed over a finite area in the complex plane, whereas the zeros of the partition 

t> ■ 

Q\ • function in the complex mass plane remain on a curve. In this paper we study the effects 

of the fermion determinant at nonzero chemical potential on the Dirac spectrum by means 
of the resolvent, G(z), of the QCD Dirac operator. The resolvent is studied both in a 
^ ■ one-dimensional U(l) model (Gibbs model) and in a random matrix model with the global 

symmetries of the QCD partition function. In both cases we find that, if the argument 
$_i ■ z of the resolvent is not equal to the mass m in the fermion determinant, the resolvent 



diverges in the thermodynamic limit. However, for z = m the resolvent in both models 
is well defined. In particular, the nature of the limit z — > m is illuminated in the Gibbs 
model. The phase structure of the random matrix model in the complex m and /i-planes 
is investigated both by a saddle point approximation and via the distribution of Yang- 
Lee zeros. Both methods are in complete agreement and lead to a well-defined chiral 
condensate and quark number density. 



March 11, 1997 



1. Introduction 



During the past decade, numerical lattice QCD simulations have provided a great deal 
of insight regarding strong interactions at finite temperature (see fl]|| f° r a review). In 
spite of steady progress numerical simulations of QCD at finite baryon density have 
proved to be remarkably difficult and their results remain inconclusive [f|,f|. The main 
reason for these difficulties can be traced back to the way the quark number chemical 
potential enters a Euclidean field theory. It results in a fermion determinant and a prob- 
ability measure that are not positive definite, making standard Monte Carlo simulations 
impossible. A possible solution to this problem is the quenched approximation, i.e., ig- 
noring the fermion determinant in generating gauge field configurations. This method has 
been quite successful in other applications of LGT (see e.g., ||). Unfortunately, results 
obtained for finite chemical potential show a critical value of n at half the pion mass 
instead of a third of the nucleon mass [|[. Two possible explanations have been offered 
in the literature: i) The problems are an artifact of relatively small lattices and would 
be cured in the continuum limit f?]||. ii) The quenched approximation is responsible for 
these unphysical results [^,10|. Regarding the latter point we mention one particularly 
interesting suggestion [11,I2|, namely that the quenched approximation is obtained as 
the limit of the number of flavors going to zero with an equal number of fermions and 
conjugate fermions, i.e., the limit of a partition function in which only the absolute value 
of the determinant enters. Recently, this suggestion has been made much more explicit 



by Stephanov within the framework of a random matrix model with the chiral and 



flavor structure of the Dirac operator (see [IS] for a review). He has shown analytically 
that the quenched Dirac spectrum is obtained in the limit as both the number of quarks 
and conjugate quarks tend to zero. Remarkably, the chiral condensate of this model (with 
the full determinant included) is independent of the number of flavors if Nf is a positive 
integer. However, since real QCD does not have such conjugate quarks, it seems unlikely 
that its Nf —>■ limit would lead to the quenched approximation. 

This leads to the more general question regarding the extent to which the quenched 
approximation can be trusted. Another familiar example of its failure is the chiral limit, 



in which the quenched theory does not reproduce the correct chiral logarithms The 
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failure of the quenched approximation is most obvious if the quark mass is less than 
the smallest nonzero Dirac eigenvalue. Then the chiral condensate in the full theory is 
zero for two or more flavors, finite for one flavor and possibly divergent in the quenched 
approximation f!3 |. The reason for the failure is that the leading contribution to the 



expectation value of an operator is a subleading contribution to the partition function. A 



similar mechanism is at work in one-dimensional U(N) lattice models p3yi5|| . There, the 
leading contribution to the fermion determinant cancels after averaging over the gauge 
group. This is no longer true if the gauge group is SU(N). Then, the chiral condensate 
in the quenched and unquenched theory are indeed equal 0,0. 

The random matrix model studied in this paper is the extension to nonzero chemical 
potential of the model introduced in [ 17| ■ The partition function is defined by 

Z = J DHP(H)det N f(m - D) , (1) 

where D is the random Dirac operator 



and W is a complex matrix distributed according to the Gaussian probability distribution 
distribution P(W) 

P(H) = exp(-nE?TrWW*) . (3) 

For nonzero chemical potential, the Euclidean QCD Dirac operator is D = ^^(d^ + iA^) + 
7o fi. The random matrix partition function is obtained from the QCD partition function 
by replacing the /i-independent terms of the matrix elements of the Dirac operator in a 
chiral basis by Gaussian random variables. In this basis, the term /i7 results in the term 
/j, times the identity in the off-diagonal blocks in (||]). We wish to emphasize that this 
random matrix model is a schematic model of the QCD partition function. 

The random matrix model ([!]) is in the class of nonhermitean random matrix models. 
Recently, a variety of such models with applications to different physical problems have 
been discussed in the literature (see for example []20| , |2l| ). 

The Euclidean Dirac operator at non-zero chemical potential is non-Hermitean, and, in 
general, its eigenvalues will be distributed in the complex plane. Indeed, this was found 



in lattice simulations by the Urbana group Q . The quenched case is defined as the theory 
which is obtained by ignoring the fermion determinant in the partition function (|l|). In 
[12 1 it was pointed out that this theory is not obtained by taking the limit Nf — * in (P, 
but rather from the partition function with det(m — D) N f replaced by |det(m — D)\ f. 
The partition function ([!]) will be studied with the help of the resolvent defined by 

G(z) = —(Tr^—) (4) 

where the average is over the probability distribution (|3|) and the fermion determinant. 
For z = m, the singularities of the resolvent will be cancelled by the fermion determinant 
in ([I]). This is the physical (unquenched) case. 

In this paper we will investigate the resolvent for z^mat nonzero chemical potential. 
Our main point is that quenching completely changes the result for G(z). In particular, 
we hope to convince the reader that, below a critical value of z, the thermodynamic limit 
of G(z) does not exist if the valence quark mass is (i.e., z) is different from the sea quark 



mass m. We shall show this for the U(l) model introduced by Gibbs | 1Q| . For the random 
matrix model ([!]) the resolvent can be obtained analytically only for z = m. In this case, 
it is well defined and displays a first-order phase transition at a nonzero critical value of 
fi. However, if valence and sea quark masses are different, numerical evidence shows that 
the thermodynamic limit of the resolvent is divergent in this model as well. 

We start with some general definitions and discuss the relations between the Dirac 
spectrum and the zeros of the partition function. Two simple models that lead to nonan- 
alyticities in G(z) are discussed in section 3. In section 4 we analyze the resolvent in the 
Gibbs model P,|j"U[. In section 5 the random matrix model is analyzed by means of 



a saddle point approximation and via the Yang-Lee zeros of the partition function (some 
of these results have been published in The quark number density is discussed in 



section 6. Numerical results for the random matrix resolvent with z ^ m are presented 
in section 7, and concluding remarks are made in section 8. 
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2. Generalities 



We will study the spectrum of the Dirac operator via the analytic properties of the 
resolvent (Q). All eigenvalues of the Dirac operator @ occur in pairs ±A. Therefore, the 
resolvent satisfies 

G(z) = -G(-z) . (5) 

If the distribution of the phases of the matrix elements of the Dirac operator over the 
ensemble is reflection symmetric with respect to the real axis, we also have the relation 

G{z*) = G*{z). (6) 

Contrary to (|5[) this relation is in general only valid after averaging. For a Hermitean 
Dirac operator, both relations are valid before averaging. 

Our central object of interest is the chiral condensate S defined by 

S = lim 5(m) , (7) 

where 

S(m) = lim — ^— d m \ogZ(m) . (8) 

The differentiation with respect to m can be carried out before or after averaging over 
the gauge field configurations. In the first case, we factorize the fermion determinant as 

det(£> - m) = Y[(m - \ k ) (9) 

k 

resulting in 

According to (|J), a nonzero chiral condensate necessarily implies that G(m) shows a 
discontinuity at m = along the real axis. 

In the second case, we factorize the partition function as 

Z(m, /i) = Y[(m - m k ) , (11) 

k 



where the rrik are the Yang-Lee zeros of the partition function. In this case 

E(m) = -J— V . (12) 

v ; N f N Y m-m k K ! 

In the phase of broken chiral symmetry, S(m) shows a discontinuity at m = 0. This 
necessarily implies that the zeros form a cut in the thermodynamic limit. 
The quark number density is defined as 

n « = WN f d » hgZ - (13) 

Here, too, we can differentiate before or after averaging over the gauge field configurations. 
In the first case, the quark number density is given by 



1 1 



n. 



q NNfff-Ck 

where the (k are the eigenvalues of 70 (D + m). They are the analogue of the eigenvalues 
of the propagator matrix introduced by Gibbs || . If, on the other hand, we first average 
over the gauge field configurations, the partition function can be factorized as 

Z(m,fi) = Y[(fi - fi k ) (15) 

k 

and results in the quark number density 

1 



n, 



~ NN f VM-Mfe 



Physically, we expect that the quark number density is zero below a certain nonzero 
critical value of \i determined by the baryon mass. This can be achieved if the eigenvalues 
(k are distributed with axial symmetry in an annulus in the complex plane. Indeed, this 
is what has been observed numerically for QCD by Gibbs. Generally, we expect p2|| 
that the zeros of the partition function are located on a one-dimensional manifold in the 
complex plane. In order to have a zero baryon density below fi c , a natural expectation is 
that the zeros are distributed homogeneously along a circle with radius fi c . 

The resolvent is analytic in the upper complex half plane for a Hermitean Dirac op- 
erator. For nonvanishing chemical potential, the eigenvalues are scattered in the com- 
plex plane. Denoting the real and imaginary parts of the eigenvalues by A r and Aj, the 



eigenvalue density is characterized by a two dimensional spectral density p(A r ,Aj). It is 
normalized according to 



J dX r d\ip(X r , Xi) = 1 . (17) 



Using the fact that 



d- z - z = tt6 2 (z) 



where the complex delta function is defined as S 2 (z) = 5(Ke(z))5(lm(z)), we find 

(19) 



p(A) = -d- z G{z) 

7T 



z=X 



One of the aims of the present paper is to understand how G(z) develops nonanalyticities 
in the complex z-plane. 

Clearly, G(z) is an analytic function of z outside the support of the spectrum of the 
Dirac operator. All singularities must be inside the support. In general, G(z) will have 
cuts in the complex plane localized within the support of the spectrum. Because of fll9|), 
the resolvent will be a function of both z and z on the support of the spectrum. 

The resolvent can be interpreted naturally in terms of a two-dimensional electrostatic 
problem ||. The electric field is given by the real and imaginary parts of the resolvent 
and the spectral density can be interpreted as the charge density. Eq. (|i~9"D is then the 
two-dimensional analogue of Gauss's law: 

1 

2n 



p(X r , Xi) = — (d x ReG(z, z) + d y lmG(z, z)) 



(20) 

(x,y)=(\ r ,\i) 



3. Simple Models 

3.1. Eigenvalues distributed homogeneously in the complex unit circle 

To illustrate the appearance of nonanalyticities in a function that shows only an explicit 
dependence on z, we consider the resolvent of eigenvalues distributed uniformly inside the 
complex unit circle. Using the electrostatic analogy we can conclude without calculation 
that 

G(z) = 0(|z|-l)- + 0(1 -1*1)* • ( 21 ) 

z 

Let us see how we can understand this result using the average over the spectral density 

G(z) = I d 2 X-^— . (22) 

J unit circle Z — A 



In polar coordinates, this integral can be written as 

G(z)= frdrf - r J -^, (23) 
Jo Jc —i {zu — r) 

where the contour integral is along the complex unit circle. The integrand has a simple 

pole at u — r/z. If \z\ < 1, we obtain a contribution only for r < \z\ which results in the 

nonanalyticity mentioned above. 

3.2. Solution of the 2x2 matrix problem 

In this section we study the resolvent for the case when 73 is a 2 x 2 matrix. The matrix 

z — D can be inverted and leads to the average resolvent 

z r°° 2 c 2v I 

G(z) = - / tdte-* / d<p- - - -— - — . (24) 

w n Jo Jo y z 2 + t 2 - fi 2 - ifit^e** + e-**) y ! 

We write the angular integral as a contour integral along the complex unit circle 

G(z) = — / tdte" 1 \ —— . 25 

w TxiJo Jc u(z 2 + t 2 - /i 2 ) -ifd(l +u 2 ) y ' 

The poles of the integrand are given by the roots u\ and u 2 of u[z 2 —t 2 + fi 2 ) —ifit(l+u 2 ) = 
0, 

. z 2 + t 2 -^ 2 ( (z 2 + t 2 -^ 2 ) 2 \ 1/2 

where u\ is the root with the plus sign. Clearly, the product of the two roots is 1. Therefore 
nonanalyticities can only appear if the roots cross the unit circle, i.e., for \u\^\ = 1. From 
the expression for the roots we can see that this happens if 

z 2 -)_ t 2 — jjL 2 

is purely imaginary, and 

2 fit 

a 2 

< 1 . (27) 



z 2 + t 2 - fl 2 



2fd 

With the aid of the real quantity t c defined below, these conditions are equivalent to 

i 

,2 j.2 2 / 2 i — 2 \ 

t = t c = fi - -(z +Z ), 

\Rez\ < fi . (28) 

At the critical point the two roots interchange resulting in the resolvent 

r°° oil r*c o 1 I 
G(z) = 2z / te- f Az6U - |Rez|) / te~* . (29) 

JO fit U\ — 112 JO fit U\ — U 2 



The first term is analytic in z\ the second term, (through t c ) depends both on z and z. 
The spectral density is then given by 



p(X r , Xi) = -d s G(z) 

71 



z=X r +i\i 



e{fi-\K\) 



1 (X 2 + X 2 )e'^-^) 



7T 



^ - xi j y + a 



(30) 



It is straightforward to repeat this calculation for the unquenched case. The condensate 
can be obtained from the partition function which, for one flavor, is given by 



Z(m)= / tdte~ l / d(j)det(m-D) 
Jo Jo 



where 

det(m — D) = m 2 + t 2 — jj 2 — ifjct^ + e" i0 ) 

The integrals are elementary resulting in the partition function 

Z(m) = n(m 2 - fi 2 + 1) . 

The corresponding chiral condensate is given by 

\ 2m 



fi 2 + l 



For one flavor the resolvent is defined by 

G N f =1 ( z ) = —- tdte' 1 \ d6-Tr det(m-D) 

v ; Z(m) Jo Jo 2 z — D ' 

which can be written as 

G N f =1 (z) = - 



fi 2 + l 



z+(m 2 -z 2 )G N f =0 (z) 



(31) 



(32) 



(33) 



(34) 



(35) 



(36) 



We see that nonanalyticities are absent for z = m. This was to be expected because 
the singularities of the resolvent are cancelled by the zeros of the determinant. A similar 
cancellation takes place in the Gibbs model to be discussed in the next section. 



4. The resolvent in the Gibbs model 



In this section we analyze the resolvent of the model proposed by Gibbs fTO ]. In this 
model the Dirac operator is defined by 



(37) 



where the indices are modulo N. Anti-periodic boundary conditions result in an extra 
minus sign for the matrix elements Df n and D G X . The partition function is obtained by 
averaging over 9, 

Z = J d9det N f {D G + m) . (38) 
The eigenvalues of D G are given by 

27ri(fe+l/2) , .„ , 27ri(fc+l/2) . a 

\ k = e K n' } +*+i* - e 4t^-^-m ? (39) 

where k = 1, • • • , N. Geometrically they fall on an ellipse in the complex plane with 
semiaxis e M — e~ M in the real direction and semiaxis e M + e~ M in the imaginary direction. 
The fermion determinant follows from the solution of a recursion relation 

det D G = e- N{i6+ ^ + e~ N ^ + X 2 (m) N + X^mf , (40) 

where 



z L z 2 



\ l , 2 (z) = -T^l + j ■ (41) 
The resolvent, defined by 

will be evaluated for Nf = and Nf = 1. (Here, Z is the partition function) 

The integral over 9 can be performed conveniently by contour integration. The phase 
of the eigenvalues and the chemical potential can be absorbed in 9 resulting in an integral 
over a circle in the complex plane with radius e M . For Nf = 1 we find 

G{z) = 9{\\ 2 {z)\-e») — -J— - (l- ^ 



X 2 (z)-X 1 (z)\ \?(m) + \$(m)J 

This is valid with the convention that |A 2 (;z)| > | Ai(a?) j . 

Note that |Ai| = \X 2 \ = 1 when z = is is purely imaginary and s < 2 and positive. For 
s > 2 our convention reads 

is i /— 



is 1 /-^ 

= 2 + 2 V ^ 1 - 

(44) 
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The condition |A 2 (z)| = e M can be rewritten as 

z = e M — for z real 

z = e M + e~ M for z imaginary . (45) 

This implies that z coincides with the modulus of the eigenvalues on the real and imaginary 
axis, respectively. 

In the quenched approximation, we obtain 

G(z) = 9(\X 2 (z)\ - e») 1 , (46) 

A 2 (2:J - Ai(zJ 

which is valid for the same convention, |A 2 (z) | > |Ai(z)|. 

We find that quenched and unquenched results agree in the thermodynamic limit if \z\ 
is larger than the modulus of any of the eigenvalues. For N — > oo and real z = r inside 
the ellipse of eigenvalues, the unquenched resolvent diverges for m < r < e M — e~ M and is 
zero for r < m. The quenched result is zero. Quenching works if m is outside the ellipse 
of eigenvalues. For imaginary z = is, quenching succeeds for s < \Jm 2 + 4 but fails in the 
region \/m? + 4 < s < e M + e~ M . For m > — e _/1 quenching works for all values of s. 

In Fig. 1 we illustrate the above quenched and unquenched results for the case m — 0, 
[i = 0.2 and N = 8. We show the resolvent for z real (left) and z imaginary (right). This 
figure clearly shows the region where quenching fails. 

It is instructive to consider the thermodynamic limit of G(z) for z — > m in the un- 
quenched case. We expect that this limit exists because the singularities of the resolvent 
are cancelled by the zeros of the determinant. Indeed, for z — > m we find 

G(z) = 9(\X 2 (z)\ - e") 1 + ^(e^ - |A 2 (*)|) ? / w S f , (47) 

A 2 (m) - Ai(ra) A 2 (m) - Ai(m) \\ 2 {m) J 

which is finite for z — m. The result for G(z) is independent of z and identical to 
jjd m logZ in the thermodynamic limit. 

5. Random matrix model 

In this section we study the partition function 

Z(m, „) = / DWP(W)det N f ( ^ ^ ^ ) , (48) 

11 



8 

Re(G) 
6 
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1.0 T , ,2.0 3.0 
Im(z) 
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0.0 0.2 0.4 0.6 0.8 1.0 0.0 
Re(z) 

Fig. 1. The real (left) and the imaginary (right) parts of the resolvent in the 
Gibbs model on the real and imaginary axis, respectively. Results are shown 
for \x = 0.2 and n = 8. The full curve represents the result for one flavor and 
the dashed curve for zero flavors. 



where W is an arbitrary complex nxn matrix and DW the Haar measure. The probability 
distribution P(W) is given by 



P(W) = exp(— n 

The fermion determinant in ( [48|) can be written as a Grassmann integral. 



(49) 



Z(m,fi) = J VWD^Vif) exp 



Nt 



k=l 



exp[-n£ 2 TWW 4 ] 



This enables us to perform the W integration 



+ v(^1 + ^l^L)] ■ (51) 
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The four-fermion terms can be written as the difference of two squares. Each square can 
be linearized by a Hubbard-Stratonovitch transformation according to 

r a 2 

exp{-AQ 2 ) ~ / daexp{-—-iQa) . (52) 

Thus, the fermionic integrals can be performed, and the partition function can be written 
as a single integral over the complex Nf x Nf matrix, a, 

Z{m,n) = J 2*7exp[-n£ 2 lW]det" ( ° + ™ fft + ro ) ■ ( 53 ) 
The condensate is given by 

(qq) = 77\rd m \ogZ^G(m) . (54) 
ZnTsf 

For n — > oo, it can be evaluated with the aid of a saddle point approximation. The saddle 
point equations are given by 

- Y?o + (a + m) ((a j +m)(a + m) - fx 2 )' 1 = , (55) 



and an equation with a and a* interchanged. The solutions are proportional to the identity 
with diagonal elements given (for S = 1) by 



a = a 

\2 ,.2 



aim + a) — fi a = m + a . (56) 



This equation was first derived p3| for the finite temperature version of the model [p3~|25||. 



(obtained from fl53|) by the substitution fi — > —iT). The condensate is then given by 

(qq) = G(m) = a . (57) 

The quark number density defined in (0) can be related to a, 

n q = (58) 

For m = it is possible to have a discontinuity in (qq) but not in n q . Notice that in the 
restored phase with o = this expression has to be treated with care. In fact, from the 
saddle point equation we obtain n q = for a — 0, and n q = —jjl for a ^ 0. In general, 
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we expect that a discontinuity in the chiral condensate is accompanied by a discontinuity 
in the quark number density 

Two solutions of ( |56|) coincide if the discriminant of the cubic equation, 

D 3 = ± (mV 2 - m\2^ - ^ - i) + (1 + /x 2 ) 3 ) , (59) 

vanishes. 

In the finite temperature version of the present model [23[|, we evaluated the average 



resolvent for z on the positive real axis and obtained the resolvent for Re(z) > through 
analytical continuation. This is justified because the resolvent is analytic in this part of 
the complex plane. The average resolvent was calculated with two methods (i) the 
super-symmetric method, and (ii) a self-consistent equation for the resolvent for n —>■ oo. 
The first method is much more intricate with regards to analyticity and convergence than 
the second. For the moment, we discuss only the second method. The self-consistent 
equation was obtained from a series expansion of the resolvent G(z) in powers of z~ l . 
This is justified if \z\ is larger than the largest eigenvalue. In our present model with a 
chemical potential, this is still true. It should also be clear that the resolvent is analytic 
in z everywhere outside of the support of the eigenvalue spectrum of the Dirac operator. 
In other words, the singularities of G(z) should be inside the support of the spectrum of 
the Dirac operator. However, for z inside the domain of eigenvalues, the series expansion 
is not justified. Indeed, the resolvent in this region is not given by the cubic equation 



(56), as was shown by Stephanov [|T^]. After replacing T — > ipL and making an appropriate 



rotation in the complex plane (the Dirac operator in @,[26[ differs by a factor i from the 
Dirac operator in the present work), the self-consistent cubic equation in p3|j2B| is seen 
to be identical to the cubic equation (|56|). 

5.1. Unquenched partition function 

In this section we analyze the partition function fl53|) for one flavor. In units where 
S = 1, it is given by 

Z(m, /i) = J dada*e- na \aa* + m{a + a*) + m 2 - /i 2 ) n . (60) 

After shifting the real part of a by m, we use polar coordinates as new integration vari- 
ables. If we notice that the angular integral represents a modified Bessel function, the 
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partition function can be written as 

Z(m, n) = 7Te~ nm / du(u - \i 2 ) n U2mn^u)e- nu . (61) 
Jo 

In the thermodynamic limit, this partition function can be evaluated by a saddle point 
approximation. Using the asymptotic form of Io(z) ~ e z /\ / 2tiz, the saddle point equation 
reads 

1 m 



u — jj? y/u 



(62) 



A phase transition takes place at the points where \Z u=Ub \ = \Z u=Ur \, with «& and u r two 
solutions of the saddle-point equation fl6"2"l) . This condition can be rewritten as 

|K - n 2 )e 2m ^- Ub \ = |(/i 2 - u r )e 2m ^- Ur \ . (63) 

The selection of the dominant saddle-points requires a detailed analysis of the partition 
function (see below). This will lead to some restrictions on the range of \x and m deter- 
mined by the zeros of the discriminant (|59|). Below we will show that (E3) determines the 
location of the zeros of the partition function. The additional limitations on the range of 
the critical values of m and fi immediately follow from such analysis. 

In general, the solution of these equations is cumbersome. However, for m^Owe find 
that u r = and Ub = 1 + fx 2 . This leads to the critical curve 



Re 



l + /i 2 + log/i 2 =0 . (64) 



For real /i, the solution is given by \i c = 0.527 . . .; for purely imaginary //, we find \x c = i. 
Moreover, from a detailed saddle-point analysis it can be shown that the partition function 
only shows a discontinuity for those values of [i that in addition to ( |64T) also satisfy < 1. 

We wish to analyze the zeros of the partition function ( |60"D in the complex m plane 
and the complex fi plane. In order to obtain the polynomial in m and /i, we expand the 
binomial in (|61~D and perform the integration over u. The result can be expressed using 
confluent hypergeometric functions 

Z(m, /,) = ^^r- E -JT- iF 1 ((n~k + 1), 1, m 2 n) . (65) 

^ k=0 

According to the Kummer identity we have 

e'^xF^in - k + 1), 1, m 2 n) = iF x (-(n - k), 1, -m 2 n) (66) 
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which is a polynomial of order n — k in —m 2 n. From the standard representation of this 
Kummer function we obtain 

I n n—k -r / 7 \ 

%ri = E E ( 7 ) • (6T) 

Although the sums in (|67|) are finite, the alternating nature of the series and the degree 
of cancellation make numerical evaluation difficult. 

For real \i the phase structure of the partition function (|6l| ) can be clarified by the 
decomposition 

Z(m, n) = Z hmken (m, //) + Z restored (m, , (68) 

where 

Z brokcn (m,/i) = 7re- nm / du(u - ^ 2 ) n I (2mn^e- nu 
Z restored (m,/i) = Tre- nm2 T du(/i 2 -u) n I (2mn^)e- nu . (69) 



o 

From a graphical representation of both sides of the saddle point equation (|62|) it is clear 
that one solution is in the interval [0,/i 2 ], which we call u r , another is in the interval 
[/i 2 , oo), which we will call Ub and another has y/u < 0. A phase transition takes place at 
the point where \Z u=Ub \ = \Z u=Ur \. For m = the thermodynamic limit of the partition 
function is then given by 

3/2 3/2 

Z(m = 0,/i) = 9(fi c -^^e-^ l+ ^ + 9^-ii c )^^ n . (70) 

'n \/n 



The partition function can be calculated in a much more effective way if we rewrite it as 
a sum of positive definite terms. For the broken part, we expand 1$ in a power series and 
change integration variables according to u — > /i 2 (l + u). The integrals can be performed 
following expansion of the binomial and result in 

z«- (m , rt = _ g g (™ 2 «)'(A)'- ■ (71) 

For the restored part, we also expand the modified Bessel function but change integration 
variables according to u — > fi 2 (l — u). The integral of each term in the expansion of 
exp(n/z 2 n) is a beta function. Collecting powers and factorials we obtain 

d (m, M ) = V (n+1) e- (m + ^ } E E „ , r [ 1 n , (m 2 n)'(A) <+a • (72) 

This expression has been used for a numerical study of the partition function. 
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5.2. Yang-Lee zeros 



In this section we evaluate the Yang-Lee zeros of the partition function ( 67p and show 
that their location is consistent with a saddle-point analysis of the partition function. 
First, we consider the partition function as a polynomial in m. Results for \x = 0, 
jj, — 0.5, and fi = 0.6 are shown in the left column of Fig. 2. We have calculated the zeros 
for different values of n, e.g., n = 48, n = 96, and n = 192. The results for n = 192 
are represented by the points in the figure. Of course, the exact location of the zeros is 
extremely sensitive to numerical round-off errors. Thus, the present results were obtained 
with the help of a multi-precision package [ 28| . Typically, we performed our computations 
with 100-500 digits accuracy. 



The zeros fall on a curve and are regularly spaced^]. It is clear from the saddle point 
analysis of the partition function that a single condition is imposed on the the complex 
variable m by the condition that the free energy of two different saddle point solutions 
coincides. This explains the fact that, in the thermodynamic limit, the zeros are located 
on a curve in the complex plane. (A similar argument has been given for the Ising model 
||22|| .) If we increase the order of the polynomial by a factor of 2, we find that half of the 
zeros are close to those of the lower-order polynomial. The other half of the zeros are 
roughly mid-way between adjacent zeros of the lower-order polynomial. This leads us to 
the conclusion that the zeros become dense and lead to a cut in the complex m plane in 
the thermodynamic limit. 



The stars in Fig. 2 represent the points at which the discriminant (|59"D of the cubic 
saddle point equation vanishes. These points coincide with the endpoints of the line of 
zeros. In the right column of Fig. 2 we show the curves in the complex m-plane across 
which there is a transition in the dominant solution of the third-order equation (^). The 
analytical result for this curve is given in (|63]). A schematic picture of the latter was also 
shown in ||29[| . 



In the left column of Fig. 3, we show the zeros of the partition function in the complex 
plane for n=192 and masses m = 0, m = 0.3, and m = 1.0. We have evaluated the zeros 



1 If the numerical accuracy is not sufficient, one typically observes that the line of zeros ends in a circle. 
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Fig. 2. The zeros of the partition function in the complex m plane (left), and 
the curves where the mean field result for the resolvent shows a discontinuity 
(right). Results are given for /j, = (upper), \x = 0.50 (middle) and fx = 0.60 
(lower) for n = 192. Zeros of the discriminant of the cubic equation ( |56| ) are 
denoted by stars. 



for n = 48 and n = 96. Here, too, all calculations were performed with 100-500 digits 
accuracy. The zeros are regularly spaced, and their density increases homogeneously with 
n. In the thermodynamic limit, we therefore expect that they join into a cut. In all 
cases, the line of zeros terminates at a zero of the discriminant (|59|), which is denoted by 
a star. In right column of Fig. 3, we show the curve across which the saddle point result 



for the chiral condensate shows a discontinuity. Below (58) it was argued that this will 



lead to a discontinuity in the quark number density as well. Notice that for m = the 
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Fig. 3. The zeros of the partition function in the complex // plane (left) and the 
location of points where the mean field result for the chiral condensate shows 
a discontinuity (right). Results are given for for m = (upper), m = 0.30 
(middle), and m = 1.0 (lower) for n = 192. The stars represent zeros of the 
discriminant of the cubic equation (|56|). Note that the scale on the x-axis of 
the lower figure is different. 



chiral condensate jumps from a finite value to zero. The analytical result follows from the 
transcendental equation fl63| ) For m = 0, chiral symmetry is broken in the region enclosed 
by this curve and is restored outside. In this case, the density of zeros approaches zero 
near fi — i. This is not surprising since, at this point, the phase transition changes from 
first to second order. 
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6. Quark number density 

In this section we study the effect of the fermion determinant and eigenvalue correlations 
on the quark number density defined in (|i~3l). In the quenched approximation the quark 
number density is given by 

nq= N Tr { tW m fI ' iW™+ii) ' (73) 
where W is distributed according to (^). For m equal to zero, the quark number density 
follows from the spectral density of the ensemble of matrices 



iW 
iW* 



(74) 



This ensemble was analyzed by Ginibre p7| . An arbitrary complex matrix can be diago- 
nalized by a similarity transformation 

W = XAX- 1 , (75) 

where A is a complex diagonal matrix. In this case, the Gaussian probability distribution 
depends both on A and X . Thus, the problem of determining the joint probability distri- 
bution of the eigenvalues is not solved by simply choosing the eigenvalues and eigenvectors 
as new integration variables. However, using a remarkable precursor of an Itzykson-Zuber 
integral, Ginibre was able to solve the problem analytically. In the thermodynamic limit, 
the eigenvalues of W (and therefore of W\ iW and iW*) are distributed uniformly in the 
complex unit circle. As we will show below, results in a nonzero quark number density 
for arbitrary small \x. In terms of the eigenvalues, the baryon number density is given by 

1 1 



1 E 



(76) 



For eigenvalues distributed uniformly inside the unit circle, the baryon number density 
can be evaluated in the same way as the chiral condensate in section 2. The result is: 

n ? = 0(l-//)/i + 0(/i-l)- . (77) 

For the unquenched case, the baryon number density in the broken phase is quite different 
with 

n q = 9(fj, c - //)(-//) +6(fJL c - fj)-. (78) 

fi 
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Finally, we wish to stress that the correlations of the eigenvalues play an important role 



in establishing the result For uncorrelated eigenvalues distributed homogeneously 

inside the complex unit circle, the unquenched partition function factorizes into 

1 



zoo 



- / d 2 X(fi 2 - |A| 2 + 2/i(A + A*)) 

7T JD 



The quark number density, which is then given by 



n„ 



2 1 ' 



(79) 



*0) 



shows some of the characteristic features of the unquenched result, including a disconti- 
nuity at a critical value of \i. In terms of the joint eigenvalue distribution of the random 
matrix ensemble (|l|), the partition function is given by 

Z (») = "X7 / IT ^A fc cZA^ n[/x 2 — |A fc | 2 -h i/x(A fc -K A%)] II |A fc — A z | 2 exp(— |A| 2 ) , (81) 

J k k k<l 

where M is a normalization constant. The correlations between the eigenvalues which are 
mainly induced by the Vandermonde determinant in (|81] leads to a partition function that 
is completely different from (Blj). (In the thermodynamic limit ( [51]) reduces to ([70D). In 
particular, the zeros of flSl|) are n-fold degenerate, whereas the zeros of (|60|) are located 
on a closed curve as shown in Fig. 3. 

It is equally straightforward to obtain the result with the determinant replaced by its 
absolute value. Then, 

1 



71 



( d 2 X(n 2 + |A| 2 + vi(A + A*)) 

JD 



(/i 2 + 



The quark number density, now given by 



n„ 



/i 2 + | ' 



(82) 



is much closer to the quenched result. 

We conclude that quenching has a major effect on the quark number density. We have 
also found that correlations between the eigenvalues play an important role in establishing 
the /^-dependence of the partition function. 
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7. Numerical results for the unquenched resolvent 

In section 4 we demonstrated that, for z less than some critical value, the resolvent of the 
unquenched Gibbs model diverges in the thermodynamic limit. We note that this model is 
special in the sense that its eigenvalues fall on a one-dimensional manifold in the complex 
plane (i.e., an ellipse), which is not the case in lattice QCD simulations 0. In this section 
we study the model ([I]), which has a genuinely two-dimensional spectral distribution. 
Unfortunately, we have not been able to treat the unquenched case analytically. However, 
it can be attacked numerically by brute force, i.e., the fermion determinant is not included 
in the measure but rather in the observable. Results for G(z) on the real and imaginary 
axes are shown in Figs. 4, 5 and 6. In Fig. 4 we show the real part of the resolvent for 
\x = 0.2 and n = 24 (left) and n = 48 (right). The mass in the determinant is equal 
to zero. The full line represents unquenched results for the resolvent. The analytical 
quenched results are represented by the dotted curves. The almost horizontal dotted 
curve shows the solution given by the cubic equation ([56|). The dotted curve through zero 
is the Stephanov solution [|l^] which is given in general by 

G(z = x + iy) = ^-^-^-x-^ . (84) 

Up to corrections of order 1/n the quenched numerical result (dashed curve) is in agree- 
ment with the Stephanov solution inside the domain of eigenvalues (\x\ < 0.0732 for 
\x = 0.2 and |x| < 0.1504 for /i = 0.3) and follows the solution of the cubic equation 
outside this domain. In this region the unquenched resolvent agrees with the quenched 
resolvent up to corrections of order 1/n. The most dramatic result of this calculation is 
that the real part of the resolvent seems to diverge for z inside the domain of eigenvalues. 
This is in complete agreement with the observation made for the Gibbs model. The slope 
near zero is strictly proportional to n. This can be understood analytically from the 
absence of eigenvalues in this region. Therefore, the resolvent can be expanded as 
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The inverse moment of the eigenvalues follows immediately from the proper unquenched 
partition function (|60"D 

W^r = d m2 \ogZ = n 2 (l + /i 2 ) (86) 
2 k x k 

resulting in 

G{z) = zn{l + fi 2 ) + 0{z 3 ) , (87) 

which is in perfect agreement with the numerical results. 

Results for the real part of the resolvent for n = 24 and fi = 0.3 are shown in Fig. 6. 
We find that the maximum increases strongly with fi. This should be contrasted with the 
same model at finite temperature where the average resolvent is flavor independent. 

The imaginary part of the resolvent is shown in Figs. 5 and 6 as a function of the 
imaginary part of z. The fluctuations of the unquenched result are much stronger in this 
case. Remarkably, the unquenched result seems to agree with the quenched result (dotted 
curve) for all real z. The area of the peak near zero is proportional to 1/n which suggests 
that this is a finite size effect. Notice that inside the domain of eigenvalues the analytical 
result (i.e., the Stephanov solution) does not depend on fi. Only the point where the two 
solutions intersect is /i-dependent. 



The results for n = 48 and fi = 0.2 were obtained by averaging over one million matrices. 
Results for fi = 0.3 and n = 24 required as many as three million matrices. The difficulties 
in simulating the partition function becomes clearer if we consider the expectation value 
of the fermion determinant which is given by the partition function for one flavor. Let us 
write out explicitly the determinant in the definition (HI) of the unquenched resolvent 



N\ \z-Dj/ {{det N f (m-D)] 

where ((•••)) stands for averaging without including the determinant in the weight. The 
normalization is just the unquenched partition function. For one flavor, it was evaluated 
in (66). In the broken phase it is given by 

3/2 

Z(m = 0, u) = ^_ e - n(1+ ' l2) (89) 
in 
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Fig. 4. Real part of the resolvent G(z) as a function of Re(z) for Im(z) = 0. 




This exponential suppression makes it prohibitively difficult to calculate the unquenched 
resolvent via a Monte Carlo method. Indeed, we find that the convergence of numerical 
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simulations becomes extremely slow for fi 2 n > 1. The average resolvent is the sum of 
terms weighted with a typical value of the determinant. Because of the normalization we 

2 

expect an enhancement factor of ~ e" M . The ratios of the maximum values of the real 
part of the resolvent in Figs. 4 and 6 are consistent with this enhancement factor. 

8. Conclusions 

We have studied the role of the fermion determinant in a random matrix model at 
nonzero chemical potential. Previously, it was shown by Stephanov that the quenched 
version of this model is given by the limit Nf — > of this partition function with the 
determinant replaced by its absolute value, and that, for example, the chiral condensate 
differs from the unquenched result. We have demonstrated that the equality of the valence 
and sea quark masses is of fundamental importance in arriving at a consistent partition 
function. A saddle point analysis was shown to be in complete agreement with an explicit 
calculation of the Yang-Lee zeros of the partition function which were found to be located 
on a curve in the complex plane. If the valence quark mass and the sea quark mass are 
different, the condensate can diverge if the mass is inside the domain of the eigenvalues. 
We have shown this analytically for the [/(l)-model of Gibbs and numerically for the 

25 



random matrix model. The reason for this divergence is that, in the chirally broken 
phase, the partition function vanishes like exp(— nfi 2 ) in the thermodynamic limit. 

While we do not expect that the random matrix model considered here leads to quan- 
titative results, we believe that some of the generic features we have observed might be 
present in QCD as well. 
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